Large Scale Structures, Symmetry, and Universality in Sandpiles 
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We introduce a sandpile model where, at each unstable site, all grains are transferred randomly 
to downstream neighbors. The model is local and conservative, but not Abelian. This does not 
appear to change the universality class for the avalanches in the self-organized critical state. It 
does, however, introduce long-range spatial correlations within the metastable states. We find large 
scale networks of occupied sites whose density vanishes in the thermodynamic limit, for d > 1. 
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One of the puzzling questions about macroscopic com- 
plex phenomena concerns the mechanisms responsible for 
the large spatially correlated structures that are often 
seen in Nature. It has been proposed that Self-Organized 
Criticality (SOC) § may be one mechanism, where the 
bursty, scale-free threshold dynamics of a slowly driven 
system is intimately linked to the emergence of long- 
range spatial (and temporal) correlations in it B. An 
obvious candidate for this picture would be the stick- 
slip dynamics of earthquakes, described by the Guten- 
berg Richter power law distribution for seismic moments, 
and faults, which form a fractal pattern in the crust of the 
earth. However, the simple sandpile, or earthquake mod- 
els do not clearly show large scale structures. Further- 
more, although many macroscopic systems show bursty 
transport phenomenology, a general feature of SOC, the 
link is not yet established because questions of robustness 
and universality are not yet resolved. Since the number 
of possible models that may be studied numerically is 
inexhaustible, it is essential to determine the symmetry 
(or other) criteria for universality [0 and robustness of 
SOC. For this purpose, we consider simple models that 
may be related to analytically solvable ones. 

Here we propose what may be the simplest sandpile 
model that gives large spatially correlated structures. For 
d > 1 the avalanches in the model have a scale- free distri- 
bution with critical coefficients in the same universality 
class as the Abelian Stochastic Directed Sandpile Model 
(A-SDM) f§-§]. There it has been proven that no spatial 
correlations exist in the steady state metastable config- 
urations §. The model we introduce is closely related 
to the A-SDM. However, a change in the rule for updat- 
ing unstable sites breaks the Abelian symmetry. (The 
Abelian symmetry refers to the fact that the order for 
updating the unstable sites has no effect on the final 
state that is reached.) This symmetry breaking intro- 
duces obvious large scale structures, consisting of net- 
works of occupied sites, within the metastable states that 
are reached in the steady state, as shown in Figs. |^ and^. 
These spatial correlations are not present when the sym- 
metry is restored. The avalanches change the network 



configuration slowly, just as earthquakes change the con- 
figuration of faults slowly. During a single or a few events 
it might appear falsely that the configuration is static or 
"pre-existing" . 

Breaking the Abelian symmetry, however, has no effect 
on the critical exponents for the avalanches, for d > 1. 
There is universality and robustness for the bursty phe- 
nomena with respect to breaking the Abelian symmetry. 
Thus two systems in the same universality class with re- 
spect to the scaling behavior of avalanches show totally 
different structures of the metastable states, with one be- 
ing completely uncorrelated and the other having chan- 
nels, or networks at large scales. 




FIG. 1. The model in d — 1. All grains from unstable sites 
in row t are thrown randomly onto neighboring sites in the 
next row t + 1. 

Consider a two dimensional square lattice as shown in 
Fig. |l|. The direction of propagation is labelled by t, with 
< t < T . The transverse direction is labelled by x, with 
periodic boundary conditions. On each site, an integer 
variable z(x,t) is assigned. The i'th grain is added to 
a randomly chosen site Xi on the top row t = 0. There 
z(xi,0) — > z(xi,0) + 1. When any site acquires a height 
greater than z c it topples, transferring all the grains at 
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that site, i.e. z(x,t) — > for z(x,t) > z c . Each grain 
from a toppling site is given equal probability to go to ei- 
ther downstream nearest neighbor, independent of where 
the other grains from the toppling site are placed. For 
each toppling event, the total number of grains are con- 
served. This is true except at the open boundary t = T 
where toppling sites simply discharge their grains out of 
the system. 

Sites are relaxed according to a parallel update until 
there are no more unstable sites, and the properties of the 
resulting avalanche are recorded. Then a new avalanche 
is initiated by adding a single grain to a randomly cho- 
sen site on the top row, t = 0. An avalanche can be 
characterized by its longitudinal extent, t c , the largest 
t row affected, its width, x c , the largest transverse dis- 
tance from the avalanche origin to any site affected by the 
avalanche, its area, a, the total number of sites affected, 
its size, s, the total number of grains thrown in toppling 
events, and the maximum number of grains thrown at 
a single site which topples, n c . The fact that z is set 
to zero at a toppling site makes the model non-Abelian. 
The corresponding rule for the A-SDM for an unstable 
site is e.g. z(x, t) — > z(x. t) — 2. 

In a recent work, Dhar jjj has shown that the stochas- 
tic Manna model @], where a fixed number of grains 
are removed from toppling sites, exhibits the Abelian 
property and is a special case of the Abelian Distributed 
Processors Model. This property was used to solve an- 
alytically for the critical state properties of the A-SDM, 
since in that case the Abelian property makes the ap- 
propriately mapped dynamics invertible ||. It is only 
necessary to realize that for stochastic models, instead 
of associating probabilities with each toppling to deter- 
mine where the grains will be thrown, we can assign to 
each site an infinite stack of independent, identically dis- 
tributed random numbers. The quenched random num- 
bers in each site's stack then determine the allocation 
of grains during each toppling event. Thus, for a given 
quenched array, and initial state of the system, the order 
of updates will not change the final configuration reached 
when the number of grains thrown from an unstable site 
is fixed. This Abelian property makes the directed model 
invertible, which leads to the product measure property 
of the metastable states, and the solvability of the A- 
SDM. However, when all grains from unstable sites are 
removed in toppling then the model is not Abelian any- 
more. 

It is straightforward to generalize the definition of our 
non-Abelian sandpile model to higher dimensions, with 
the number of directions transverse to the direction of 
propagation being d. The threshold z c can be chosen ei- 
ther to scale with dimension as z c = 2d — 1 (for d > 1) 
or it can remain constant at z c = 1. The same behaviour 
and scaling exponents are recovered under both condi- 
tions. Below, unless stated otherwise, we refer to the 
model with z c = 1. 



First, we discuss the case d = 1. From the dynamical 
rules it is clear that the avalanches must, themselves, be 
essentially compact. Thus each avalanche sweeps out ar- 
eas of the lattice leaving empty sites. At the edges of the 
avalanche, sites occupied with grains may remain. Thus 
in the stationary state the structure of the sandpile will 
consist roughly of empty areas bounded by wandering 
paths of occupied sites, which can branch and recombinc. 
At the top of the sandpile, where the grains are added, 
the network of grains is dense, but pushing into the sand- 
pile it becomes coarser and coarser. This coarsening re- 
flects the fact that avalanches that reach further into the 
system are bigger and wider and thus leave traces at their 
edges that are further apart. A steady-state sandpile con- 
figuration is shown in Fig. |[ 

In fact the average density of sites occupied with grains 
scales with distance from the top of the pile as p(t) ~ t~ a , 
with a = 0.45 ± 0.02. In spite of the vanishing density in 
the thermodynamic limit, this network of grains is essen- 
tial for maintaining the steady state of SOC, providing 
a drainage outlet for grains to be transported from the 
top to the bottom of the system. The situation for the 
A-SDM is completely different. In that case, the density 
of occupied sites is 1/2 for d = 1, and the occupation 
numbers for sites are completely uncorrelated, being de- 
scribed by a product measure in the steady state || . 




FIG. 2. A steady-state configuration of our d = 1 sandpile 
model showing sites occupied with grains forming a network 
that can transport grains from one end of the sandpile to the 
other. The shaded area indicates the sites which toppled in 
the preceding avalanche. 

Despite these vast differences, the distribution of 
avalanche sizes and durations in the steady state of SOC 
exhibits finite size scaling with the same critical expo- 
nents. This is demonstrated in Fig. |^ where the size 
distributions of avalanches, for both the A-SDM and 
our non-Abelian model, are presented. All the criti- 
cal exponents characterizing avalanches in the A-SDM 
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have been determined analytically ]5|,^| and numerically 
P,|6|JT^]. Avalanches correspond to an absorbing state 
phase transition G3]. They are described by a vari- 
ant of the Edwards- Wilkinson jflj stochastic interface 
equation where the noise amplitude is a threshold func- 
tion of the height (occupancy in the sandpile model) || . 
Using these analytically determined values for critical 
exponents, good data collapse is also obtained for our 
non-Abelian model, where no analytic solution exists at 
present. Thus, within our numerical accuracy, the criti- 
cal exponents for the avalanches are the same in the two 
cases. It appears as though the non-Abelain sandpile 
has organized large scale structures in such a way as to 
maintain the universality class, governed by a stochastic 
continuum equation ||, for the avalanches. 
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FIG. 3. Finite size scaling of A-SDM and our sandpile 
model using r s = 10/7, D = 7/4. 

A configuration in the steady state of our sandpile 
model in two dimensions is shown in Fig. ^. We observe 
a type of domain tube structure with walls separating the 
different domains. The tube domains get larger as they 
go into the system. Again, a numerical analysis of the 
avalanche size distribution using finite size scaling gives 
critical exponents t s = 3/2 and D = 2, the same values 
as determined analytically for the A-SDM |]||. Simi- 
larly for the distribution of avalanche times, t, we find 
finite size scaling with exponents r t — D = 2. There is 
no such tube structure, though, in the A-SDM. 

The zero dimensional model is a chain of sites. Since 
each site has only one downstream nearest neighbor, the 
dynamical rules of the model must be specified in a 
slightly different way. We allow grains to be distributed 
to both the nearest and next-nearest neighbours down 
the chain, and consider two ways in which the relaxation 
of critical sites can be ordered. With a parallel update 
rule all critical sites are relaxed simultaneously. In this 
case, the time in terms of the parallel update at which 



a site can become critical is not equal to its row number 
and sites may topple many times during an avalanche. 
(This does not occur for d > 0.) We can also define a 
single-site update rule in which the top-most unstable 
site is relaxed each time step. Multiple toppings cannot 
occur in that case. These two cases lead to different sets 
of critical exponents for the avalanches in d = 0. 
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FIG. 4. A steady state configuration of the sandpile model 
in d — 2 dimensions. Note that the axes are not scaled pro- 
portionally and lattice sites with t < 500 have been removed 
for clarity. 

The parallel update dynamics yields the same 
avalanche exponents, e.g. D = 3/2, t s = 4/3, as the 
Abelian model that was studied by Kloster et al |(| (see 
Fig. ||). We found good data collapse for system sizes 
ranging from T ~ 10 3 to T ~ 3 x 10 4 , for both the size, s, 
and time extent, t, of the avalanches. However the aver- 
age occupancy of sites in the steady state does not decay 
to zero as the distance from the top site, where grains 
are added, increases, as in our model in higher dimen- 
sions. Instead it is constant p(t) = j, apart from small t 
where the average occupancy adjusts exponentially from 
the initial value of \. Occupied sites do not appear to be 
spatially correlated. Thus the behavior of the model is 
very similar to its Abelian relative || and appears to be 
in the same universality class. 

In the parallel update dynamics, as more than one site 
can topple at each time step and the location of the top- 
plings is allowed to vary, the avalanches themselves have 
structure. On average the active front moves through the 
system with velocity 1.5 (i.e. advances 3 lattice spaces in 
2 parallel updates, on average). Around this average the 
active sites are split into a series of smaller fronts which 
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spread, branch and recombine (see inset in Fig. ||) 

Since our model is not Abelian, a change in the rule 
for the order in which unstable sites are updated can 
have a pronounced effect. We tried a variety of update 
rules in higher dimensions d > 0, but all the versions we 
tried appeared to give the same universality class for the 
avalanches, although the structure of metastable states 
did change drastically. 

However, in d — the critical behaviour is less ro- 
bust. A finite size scaling analysis of the avalanche size 
distribution using the single site update rule described 
above does not collapse with the same exponents as the 
d = A-SDM. A reasonable data collapse can be pro- 
duced by changing the cutoff dimension to D — 1.1 and 
keeping t s = 4/3, but it is not completely convincing. 
A multifractal data collapse M did not yield noticeably 
better results. The multifractal collapse was performed 
with parameters sq = 0.5, Iq = 0.2. For the single site 
model, the density of occupied sites does decay going into 
the system from the top where the grains are added. It 
behaves approximately as p(s) ~ 1/t. 
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FIG. 5. Finite size scaling of the distribution of avalanche 
sizes for the A-SDM and our model with parallel update in 
d = 0, using t s = 4/3 and D — 3/2. Inset shows the position 
of active sites away from the average at each time step. 
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